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THE  METHOD  OF  LINES  FOR  NUMERICAL  SOLUTION  OF 
PARTIAL  DIFFERENTIAL  EQUATIONS 

ABSTRACT 

In  the  method  of  lines  for  solving  certain  kinds  of  boundary  value 
problems  in  rectangular  or  trapezoidal  regions  one  of  the  variables,  say 
y,  is  discretized  while  the  other  variable  x  is  left  continuous.  When 
suitable  finite  difference  approximations  are  substituted  for  the  partial 
derivatives  with  respect  to  y  the  differential  equation  is  changed  into 
a  simultaneous  system  of  ordinary  differential  equations  in  the  variable 
x.  Ihe  method  used  very  little  in  the  USA  is  used  extensively  in  the 
Soviet  Union  and  nearly  all  the  literature  on  this  subject  is  in  Russian. 
The  method  has  been  tried  in  BRL  and  it  seems  to  be  a  very  useful  one. 
This  report  does  not  pretend  to  be  a  monograph  on  the  subject.  It 
intends  to  be  a  practical  guide  to  computations. 
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THE  PRINCIPLE  OF  THE  METHOD  OF  LINER 

We  shall  explain  the  method  of  lines  for  the  following  differential 
equation  of  the  second  order  in  two  variables  which  is  to  be  integrated 
in  the  rectangular  region. 

R;  a  £  x  £  (3;  yQ  £  y  £  yQ  +  L 

with  boundary  C. 

Our  boundary  value  problem  is 


au 

♦  bu  +  cu  + 

du 

+  eu 

+  f  u  ■  g  in  R 

(1) 

XX 

xy  yy 

X 

y 

u(x, 

y0)  -  <J0to; 

u(x, 

l)  -  q-,00 

(2) 

on  C 

u(a, 

y)  -  p0(y); 

u(S, 

y)  - 

px(y) 

(3) 

where  a,  b,  c,  d,  e,  f,  g  are  functions  of  x  and  y,  and  q^^  and  are 
prescribed  functions  of  x  and  y  and  all  of  these  functions  are  continuous 

To  solve  the  above  boundary  value  problem  by  the  method  of  lines  we 
shall  use  the  following  procedure: 

Subdivide  the  interval  L  ■  y  -  y^  into  n  +  1  equal  subintervals 

n+l  o 

of  width  h  ■  L/(n  +  l),  then  draw  n  lines  parallel  to  the  x  axis 
y  .  yk  -  y0  +  Ishj  k  .  1,  2,  n. 

which  form  the  grid  shown  in  Figure  1. 
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We  assume  that  both  the  first  and  the  second  order  partial  derivatives 
are  continuous  in  x  and  y.  Then  we  substitute  in  Equation  (l) 

y  -  yk;  (*  -  2,  n) 

and  replace  the  partial  derivatives  with  respect  to  y  by  the  central 
differences 

Uy(x,  yk)  ~  (2»)'1[uk+l(x)  -  Uk-l(x)] 

Uyy^X)  y]j)  ~  G1)  Pk+l^X^  "  2Uk^X^  +  Uk-l^j 
v (x-  yk>  ~  (2hrl[uiti(x)  -  ^-I(x)] 

where 

Uj^(x)  «^J(x,  yk);  U£(x)  -  (d/dx)(u(x,yk)),  and  Uk(x)  is  an 
approxlaation  of  u(x,  yk)  on  the  line  y  ■  yk< 

When  we  perform  these  substitutions  we  obtain  a  system  of  n  simul¬ 
taneous  differential  equations  of  the  second  order  which  approximates 
the  system  Equations  0-)  and  (2); 

«k°k +  -  «U>  ♦  •  2Uk  *  Uk-1> ♦  Vk  w 

♦  (2h)’aek(B^sfl  -  Uk^)  +  f^  ■  g^i  k  •  1,  2,  n. 

The  boundary  conditions  become 


U0(*)  -  l^x); 

Vl(x)  *  4l(x) 

(5) 

Va)  -  r0(jrk)j 

i>k(P)  -  px(yk)  • 

(6) 

Nov  Equations  (5)  are  no  longer  considered  to  be  boundary  conditions. 
They  determine  certain  ten*  in  the  equation  for  k  ■  1  and  for  k  •  n 
belonging  to  system  (4).  The  Equations  (6)  are  the  2n  boundary  conditions 
for  the  n  second  order  differential  Equations  (4). 
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The  simultaneous  system  of  ordinary  differential  Equations  (4)  and 
(5)  together  with  the  2n  boundary  conditions  Equation  (6)  approximate  the 
boundary  value  problem  Equations  (l),  (2),  and  (3).  The  general  solution 
of  Equations  (4)  and  (5)  depends  linearly  on  the  2n  arbitrary  constants 
of  integration  which  are  determined  from  the  2n  boundary  conditions 
Equation  (6). 

The  convergence  of  the  approximating  system  Equations  (4),  (5),  (6) 

to  the  original  system  Equations  (l),  (2)  and  (3)  when  h  approaches  zero 

under  certain  restrictions  on  the  coefficients  and  on  the  boundary 

1* 

conditons  has  been  proven  by  various  Soviet  mathematicians  . 


THE  LAPLACE  EQUATION 

Consider  the  boundary  value  problem  Equations  (l),  (2)  and  (3)  when 
the  left  member  of  Equation  (l)  is  the  Laplacian 

^  "  uxx  +  Uyy  *  ‘  (**) 

In  this  case  the  approximating  system  of  Equation  (4)  takes  the  form 

U”  +  h"2(l\cfl  -  2Uk  ♦  Uk-1)  -  gk;  k  -  1,  2,  n.  (4A) 
The  Equations  (5A)  and  (6A)  would  be  the  same  as  Equations  (5)  and  (6). 

A  HIGHER  ORDER  OF  APPROXIMATION  FOR  THE  LAPLACE  (OR  POISSON)  EQUATION 

The  order  of  approximation  for  the  system  Equations  (4)  or  (4A)  is 
0(h2).  For  the  Laplace  equation  we  car.  derive  an  approximating  system 
of  the  order  O(h^).  To  obtain  it  we  expand  u^  and  tn  Taylor 
Series  about  yk,  keeping  the  fourth  order  terms,  and  after  eliminating 
the  fourth  partial  derivative  we  get 

0/6)11-  ♦  (1/12HU-J  ♦  U-. )  ♦  h-2^  -  2Uk  *  Uk.x)  m 

•  (5/6)gk  ♦  (l/lsHg^  ♦  gj^)  . 


Suptreoript  numbore  denote  re/enenoee  uhioh  may  b§  found  on  pay*  30. 
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THE  CLOSED  SOLUTION 


A  * 


For  the  Laplace  equation  and  when  the  prescribed  values  of  u  on  the  • 
lines  y  «  yQ  and  y  ■  y  ^  are  zero,  we  can  obtain  a  simple  closed  solution 
for  each  line. 

Consider  the  boundary  value  problem  Equations  (l),  (2),  and  (3), 
when  qi  ■  0,  which  is  approximated  by 

U£-  +  -  2Uk  ♦  0k  l)  -  0;  k  -  1,  2,  a  (W) 

u0M  -  Untl(x)  -  0  ,  (5A) 

u(a,  yk)  -  p0(yk);  u(P>  yk)  -  px(yk);  k  *  1,  2,  n.  (6a) 

Applying  the  separation  of  variables  we  assume  the  following  form 

Uk(x)  -  q(k)v(x) 

and  substitute  it  into  the  above  homogeneous  equation.  This  yields  the 
following  equation: 

q(k)v"(x)  h^v(x)Jq(k  ♦  l)  -  2q(k)  ♦  q(k  -  l)j  -  0 

q(0)  «  q(n  ♦  l)  -  0 
or 

V"(x)/v(x)  •  jq(k  ♦  1)  -  2q(k)  ♦  q(k  -  l)]/  -  h2q(k)  -  62  *  constant. 

Ib  find  q  we  must  solve  the  homogeneous  difference  equation 

q(k  -  1)  -  [2  -  hVjq(k)  ♦  q(k  -  1)  -  0 

with  the  boundary  conditions 

q(0)  ■  q(n  ♦  1)  ■  0  . 

The  general  solution  of  this  difference  equation  has  the  form 

l(k)  -  ♦  Cgkg 
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where  and  Cg  are  arbitrary  constants  and  X^  and  Xg  are  the  roots  of 
the  characteristic  equation 

X2  -  [a  -  h2B-]x  *1.0. 

From  the  boundary  conditions  we  have 

q(0)  »  +  Cg  ■  0,  hence  Cg  ■  - 

q(k  +  1)  -  C1(X1n+1  -  Xg1*1)  -  0,  hence  (X1/Xg)n+1  -  1 
and 

(X^Xg)  -  exp  (2fT  is/(n  ♦  l))  . 

From  the  characteristic  equation  we  have  that 

X1X2  "  1 

consequently 

Xx  -  exp(n  is/(n  +  1))  ; 

Xg  -  exp(  -  n  is/(n  ♦  1))  , 

8  -  (y8  -  y0)A  -  1,  2,  ..  n  . 

From  the  characteristic  equation  we  have  also  that 

2  2 

X1  *  X2  *  2  *  h  6 

consequently 

2  -  h262  •  exp(ti  is/(n  ♦  1))  ♦  exp(-  n  is/(n  ♦  1))  «  2  cos(tts/(n  4  1)) 
h262  »  2  -  2  cos  (ns/ (n  ♦  l))  ■  l»  •in2(n(yi  -  yQ)/2 L 

q§(k)  •  c[cxp(n  Uk/(n  ♦  1))  -  exp(-n  i*k/(»  ♦  *))]  ■  C  sin(t»s(yk  -  y0)/L). 


U 


Then  taking 

v"(x)  -  82  v(x)  «  0 

we  obtain 

i 

vg(x)  -  Csexp(8sz)  +  Dge3q>(  -  8gx)  . 

Thus,  we  have  a  set  of  linearly  independent  solutions 

ukjS(*)  "  [csexP(V>  +  Dsexp^’  *s*)]*M*s(yk  -  yo)fa)i  s  =  1,  2,  ...  n 

and  the  general  solution  is 

VX)  =  I  [csexp(8sx)  +  Dgexp(-  5sx)lsin(TTs(>k  -  yQ)/L) 
s«l 

where  C  and  D  are  arbitrary  constants, 
s  s 

In  a  similar  way  it  can  be  shown  that  the  solution  of  the  homogeneous 
syrtem  corresponding  to  the  higher  order  approximation  for  the  Laplace 
Equation  '4B)  is 

Uk(x)  =  £  [c^exp(B^x)  +  D^exp(-  6^x)]£in(ns(yk  -  yo)/L) 
s»l 

where 

6^2  -  21+  sin2(n/2L)(yg  -  yQ)/h2(5  +  cos(TT/L)(yg  -  yQ)) 

and  Cg  and  Dg  are  arbitrary  constants. 

Having  the  general  solution  of  homogeneous  system  we  may  be  able 
in  many  concrete  cases  to  find  the  particular  integral  corresponding  to 
the  given  right  member,  g(x,  y). 

For  example  if  g  is  a  constant  or  a  function  of  y  only,  then 

Au  «  g(y)  (1) 

and 

S +  ♦  Vi>  *  v  <«k  *  •  (**) 
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Assume  that  the  particular  integral  on  the  k-th  line  is 


(a  constant). 


Substituting  it  in  Equation  (4A)  we  obtain  the  linear  system  of  equations 


Vi  -  “k  +  Vi  ■Ah'2 

from  which  the  values  of  A^  are  easily  determined. 

When  g  is  a  n-th  degree  polynomial  in  x,  we  can  assume  the  solution 
in  the  form  of  another  n-th  degree  polynomial  whose  coefficients  are 
determined  by  substituting  it  in  Equation  (4A).  Let  for  example 

2 

g  a  mQX  +  DLjX  +  mg  . 

Assume  the  solution  to  be 

uk  *  v2  +  V  *  °k  • 

Substituting  it  in  Equation  (4A)  we  obtain 

A +  (Vi  -  A ♦  Vi>*2  +  <Vi  -  A +  W* 


*  Vi  -  2ck  +  °k-i  -  V2  +  V"2  • 

Comparing  coefficients  of  the  same  powers  of  x  we  have 

Vi  -  2VA-i  *  "o 

Vi  -  A  *  Bk-i  ‘  “i 
Vi  -  A  +  ck-i  + A’A 


from  which  the  values  of  V  V 


C,  can  be  determined, 
k 
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NUMERICAL  EXAMPLE  1 


Solve  the  following  boundary  value,  problem: 


Au  ■  -  1  In  the  region  R  which  is  the  rectangle  (y|  £  j|$  |x|  £  (l) 

u(|-  ,  y)  ■  u(-  I*  ,  y)  «  0  on  the  boundary  C;  (2) 

u(x,  |-)  »  u(x,  -  i)  -  0  .  (3) 


Applying  the  method  of  lines  we  shall  use  the  three  lines 
7i  ■  -  5s  y2  ■  0;  y3  -  (h  »  jp  n  -  3;  L  =  l)  . 

We  shall  compute  U^x)  on  these  lines  using  the  approximating  system 
of  Equations  (4A),  which  in  our  case  is 

U^(x)  +  16[U2(x)  -  2U1(x)]  =  -  1;  UQ(x)  =  U4(x)  =  0. 

tf^(x)  +  16[U3(x).  -  2U2(x)  +  Ur(x)]  «  -  1  (la) 


U^(x)  +  16[U2(x)  -  2U3(x)]  =  -  1  , 
with  the  boundary  conditions 


°i<*  f).VV?)  ■  08  1' 1.2-  3. 


Figure  2 
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THE  PARTICULAR  INTEGRALS 


Ve  shall  assume  that  the  particular  Integrals  of  our  system  are 
the  constants 

°i  -  Ai  • 

Substituting  the  above  in  the  system  Equation  (la)  we  obtain 
A1  =  ■  3/32;  A2  *  l/8 

THE  GENERAL  SOLUTIONS 

Since  the  prescribed  values  of  u  on  the  lines  y  *  -  —  and  y  =  are 
zero  we  can  use  the  closed  solution.  Adding  the  complementary  functions 
to  the  particular  integrals  we  can  write  the  general  solutions  on  each 
line  as 

U^x)  =  sin  jJ-  [C^exp(djX)  +  ^expC-d^x)]  +  sin  [cgexp(d2x)  +  Dgexp(-dgx)] 
+  sin(3n/4)[c3exp(d3x)  +  D3exp(-d3x)]  +  3/32  ** 

=  (/^Xc^exp^x)  +  D1exp(  -d-jX)  +  C^exp^x)  +  ^expf^x)] 

+  C2exp(d2x)  +  Dgexp(-dgx)  +  3/32  ; 

U2(x)  =  sin  |-tt [c^exp(d-jx)  +  D1exp(-d1x)]  -  sinTr[c2exp(d2x)  +  D2exp(-d2x)] 

+  sin(3n/2)[C3exp(d3x)  +  D^expC-d^)]  +  1/& 

»  C1exp(d1x)  +  D^expC-d^x)  -  C^xp^x)  -  D^expC-d^x)  +  l/8 
U3(x)  »  sin(3n/4)[ciexp(d1x)  +  D^xpC^x)]  +  sin(3TT/2)[C2exp(dgx) 

+  Dgexp(-dgx)]  +  sin(9rrA)Cc3exp(d3x)  +  D3exp( -^x)]  +  3/32 
-  (/S'/2)[ciexp(d1x)  +  D1exp(-d1x)]  -  [c2exp(d2x)  +  Dgexp(-dgx)] 

+  (y?/2)[C3exp(d3x)  +  D3exp( -djx)]  +  3/32 
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where 


d2  a  64  sin2(n/8)j  d2  ■  64  sin2(n/4);  d2  *  64  sin2(3rr/8)  . 

Since  the  region  R  and  the  boundary  conditions  are  symmetrical  with 
respect  to  the  y  axis  that  is  U^(x)  a  U^(-x),  we  have 

°i  ■  Di 

and  using  the  boundary  conditions  we  obtain  the  system  of  algebraic 
equations  which  determine 

(/§  cosh(d1/2))C1  *  (2  cosh(d2/2))C2  +  (/?  cosh(d3/2))C3  -  -  3/32 
(2  cosh(d^/2) )c^  -  (2  cosh(d3/2))C3  =  -  l/8 
[/?  cosh(d1/2))C1  -  (2  cosh(d2/2))C2  +  (/2  cosh(d3/2))C3  -  -  3/32 
Hence 

Cx  =  -  (3/2  +  4)sech(d1/2)/l28;  C3  a  -  (3/?  -  4)sech(d3/2)/l28;  C2  -  0. 
Thus  finally 

Ux(x)  a  /5cx  cosh(dxx)  +  /SC3  cosh(d3x)  +  3/32 
U2(x)  *  2CX  cosh(d1x)  -  2C3  cosh(d3x)  +  l/8 
U3(x)  a  /Scx  coshCd^x)  +  /5c3  cosh(d3x)  +  3/32 

NUMERICAL  EXAMPLE  2 

Solve  the  boundary  value  problem 

Au  ■  0  in  R  (l) 

u(x,  0)  •  u(x,8)  -  100  x  (12  -  x)  (2) 

on  C 

u(0,  y)  «  u(l2,  y)  ■  100  y  (8  -  y)  (3) 

R;  0  *  x  *  12;  0  *  y  *  8 
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y 


Applying  the  method  of  lines  we  shall  use  the  three  lines 

yl  =  Vq  ~  ^  y3  =  (h  =  2,  n  •  3,  L  »  8). 

We  shall  compute  U^(x)  on  these  lines  using  the  approximating 
system  Equation  (4b)  which  in  our  case  is 

(5/6)Uj  +  (1/12) (U1^  +  U"^)  ♦  h_2(Uk+1  -  2Uk  +  U^)  -  0  (4B) 

Uq(x)  «=  U^x)  -  100  x  (12  -  x)  (5B) 

Uk(0)  «  Uk(l2)  .  100  yk(8  -  yk);  k  »  1,  2,  3.  (6B) 

We  shall  rewrite  the  system  combining  Equations  (4B)  and  (5B) 

(5/6)U"  +  (l/l2)(u£  -  200)  ♦  (lA)(U2  -  2UX  +  100  x  (12  -  x))  «  0 

(5/6)1^  +  (l/l2)(U*2  +  Up  *  (lA)(U3  -  2Ug  +  Up  -  0  (kB) 

(5/6)tfJ  ♦  (l/l2)(-200  +  u£)  +  (lA)U00  x  (12  -  x)  -  2Uj  ♦  Ug)  -  0 


yx  ■  2j 

Ux(0)  -  1200; 

UX(12)  -  1200 

y2  * 

U2(0)  -  1600; 

Ug(l2)  -  1600 

(SB) 

y3  *  6> 

U3(0)  -  1200; 

U3(12)  -  1200 . 
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For  digital  computations  the  second  order  system  must  be  reduced  to 
an  Initial  value  first  order  system.  To  achieve  this  let  us  set 

U{-V1  ;  (I) 

-  V2  ;  (II) 

-  V3  ;  (III) 

Substituting  Equation  (i),  (II),  (III)  in  Equation  (4B)  we  obtain 

(5/6)V|  +  (l/l2)(V^  -  200)  ♦  (lA)(Ug  -  2U^  +  100  x  (12  -  x))  -  0  (IV) 

(5/6)V*  +  (l/l2)(V*  +  V’)  +  (l/4)(U3  -  2U2  +  Ux)  -  0  (V) 

(5/6)V^  +  (l/l2)( -200  +  V2)  +  (l/4)(l00  x  (12  -  x)  -  2U3  +  Ug)  -  0  (VI) 

or 

+  (l/lO)V£  +  (3/lO)Ug  -  (6/10)1^  -  20  -  30  x  (12  -  x)  (I V) 

V*  +  (i/io)(v^  +  v*)  ♦  (3/io)(u3  -  2U2  +  ux)  -  0  (V) 

V'3  ♦  (l/lO)V£  -  (6/lO)U3  +  (3/lO)U2  -  20  -  30  x  (12  -  x)  (VI) 

The  Equations  (I)  through  (VI)  form  the  system  of  six  first  order 

equations,  two  of  them  nonhomogeneous.  The  conditions  Equation  (6B), 

however,  are  not  all  initial  and  we  have  to  arrange  for  that.  The 
solution  on  any  line  k  will  be  the  linear  combination  of  independent 
solutions  of  the  homogeneous  sytem  u£  plus  the  particular  Integral  u£ 

Uk(x)  -  Y  C^(x)  ♦  U*(x)  k  -  1,  2,  3,  (VII) 

A 

where  C1  are  constants  to  be  determined  from  the  boundary  conditions. 

The  independent  solutions  of  the  homogeneous  system  Equations  (I) 
through  (VI),  where  the  right  members  in  the  Equations  (IV)  and  (VI)  are 
replaced  by  zeros,  are  obtained  from  the  following  set  of  Initial  condi¬ 
tions  at  x  ■  0: 
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TABLE  I 

INITIAL  CONDITIONS,  x  -  0 


Symbols  of  the 
Ind.  Solutions 


V°>  V°>  u3(°5  V0)  y°)  v°> 


1st  set  the  initial 
conditions 

2nd  set  the  initial 
conditions 

3rd  set  the  initial 
conditions 

4th  set  the  initial 
conditions 

5th  set  the  initial 
conditions 

6th  set  the  initial 
conditions 


The  particular  integrals  are  obtained  from  the  non -homogeneous 
system  Equations  (I)  through  (VI)  with  initial  conditions  all  zero  (as 
shown  in  the  above  table).  The  constants  C1  are  determined  from  the 
linear  system  arising  from  substituting  in  Equation  (VII)  the  boundary 
conditions  at  x  ■  0  and  at  x  ■  12  ( 


l  ciu^(o)  .  u  (0)  -  uj(o) 


l  'V12*  ■ 


1^(12)  -  Ug(l2)j  k  .  1,  2,  3. 


Vs  shall  explain  Table  1  on  ths  example*  Ths  first  lifts  indtoatss  that 
ths  initial  conditions  art  VJO)  *  1;  VJO)  -  IL(0)  -  VJO)  «  VJO)  - 
VJO )  »  0,  Ths  solutions  resulting  from  thsss  Initial  conditions  art 


and  Vi 


CURVILINEAR  BOUNDARIES 


If  the  region  of  integration  R  is  of  a  shape  of  a  curvilinear 
trapezoid  shown  in  Figure  2  which  is  hounded  by  the  lines 

y  -  y0;  «n4  y  -  yn+1 

and  by  the  curves 

x  *  a(y)  and  x  »  3(y);  y  *  y  £  y  , 

o  n-i 

then  the  procedure  of  the  method  of  lines  remains  essentially  the  same 
as  for  a  rectangular  region.  The  proof  of  convergence,  however,  requires 
that  the  third  partial  derivative  with  respect  to  y  be  continuous.  The 
curvilinear  boundaries  will  be  explained  in  the  following  example. 


EXAMPLE  3 

Solve  the  boundary  value  problem 


Au>  0  in  R 

(1) 

u(x,  0)  ■  u(x,  k)  »  0 

(2) 

on  the  curve  a  u  »  ^(x,  y)  »  x  +  y 

(3) 

on  the  curve  3  u  •  (x>  y)  ■  x  -  y 

where  R  is  a(y)  *  x  S,  3(y);  0  £  y  *  k 


x  -  a(y)  •  g-  y2  +  1;  x  -  3(y)  «  ^  ey  ♦  3 


y  • y it 


Figure  k 
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1 


2 


Applying  the  method  of  lines  we  shall  use  the  three  lines 
yl  “  1;  y2  “  2;  y3  -  3(h  -  1,  n  ■  3*  L  -  4)  . 

We  shall  compute  Uk(x)  on  these  lines  using  the  approximating 
system  of  Equations  (4A)  which  in  our  case  is 

u”  +  U2  -  2UX  +  Uq  «  0;  Uq  -  U4  -  0; 

U*2  +  U3  -  2U2  +  Ux  =  0 


tfj  ♦  U4  -  2U3  +  U2  -  w 


On  the  curve  a; 
On  the  curve  0; 


1  2 


W  *xia+yiI-ffyk  +  >,k  +  1 

"  Xk2 


k-r1'k-'k*3 


where  x^  is  the  abscissa  of  the  point  of  intersection  of  the  line  y  «  y^ 
and  the  curve  x  ®  a(y),  and  x^  is  the  abscissa  of  the  point  of  inter¬ 
section  of  the  line  y^  and  the  curve  x  ■  0(y). 

Like  in  the  previous  examples  we  obtain  six  independent  solutions 
from  the  assumed  initial  values  shown  in  Table  I  and  form  the  general 
solutions 

6 

Uk(x)  -  £  C*u£(x) 

L*1 

We  determine  the  constants  C1  from  the  linear  system  arising  from 

substituting  th*  bounuary  conditions 


V\i>  •  l  cluk(xki) 


Uk(*k2>  *  i  C‘Uk(Xk2)!  k  •  l>  2-  3- 
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MACHINE  COMPUTATION 


In  order  to  compare  the  method  of  lines  with  the  conventional  grid 
method.  Numerical  Example  2  has  been  programmed  using  the  two  methods. 

The  programming,  with  notation  included,  for  the  method  of  lines  is 
given.  Figure  5  is  a  flow  chart  showing  computer  operations. 

CONCLUSIONS 

Hie  method  of  lines  and  the  conventional  grid  methods  have  been 
compared  on  two  high-speed  digital  computers  at  Ballistic  Research 
Laboratories,  Computing  Laboratory,  Aberdeen  Proving  Ground,  Maryland, 
with  respect  to  run  time,  computer  limitations,  and  one  known  solution, 
U(U,  4)  »  2428.  Let  "H"  be  the  step  size  and  MN"  be  the  number  of  points. 
The  comparisons  follow: 

First  Method  -  Conventional  Grid  Method 

A.  H  «  2  N  «  15  15  x  15  Matrix 

OREVAC  -  Run  Time  5  min. 

No  Limitations 

BRLESC  -  Run  Time  1  min. 

No  Limitations 

U(4,  4)  -  2419.53919 

B.  H  •  1  N  -  77  77  x  77  Matrix 

ORIVAC  -  Memory  too  small 

BRLESC  •  Run  Time  1  min 

No  Limitations 

U(4,  4}  -  2420.69488 

C.  H  -  .5  N  ■  308  308  x  308  Matrix 

Memory  too  small  on  both  computers 
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nun I  ,  _  i.t  I  aouNMRY  coromcws 


Y'l-Yk  Y’U.V’l 
1*2-15  T*5-V*2 
Y*  3-T6  Y*6-V*  3 


Second  Method  -  Method  of  Lines 


A.  Is  2  N  ■  3  AX  *  .1  6x6  Matrix 

ORDVAC  -  Run  Time  5  lain* 

Limitations,  smaller  M's  consume  too  much 
run  time 

BRLESC  -  Run  Time  1  min. 

No  Limitations 

U(4,  4)  2420.4529 

B.  H  =  1  N  =  7  AX  =  .1  14  x  14  Matrix 

ORDVAC  -  Run  Time  10  min. 

Limitations,  same  as  A. 

BRLESC  -  Run  Time  1  min. 

No  Limitations 

u(4,  4)  2420.7435 

Tne  method  of  lines  needs  approximately  ten  times  less  storage  than 
the  conventional  finite  difference  methods.  In  some  cases  it  may  be 
faster  and  more  accurate.  Another  advantage  of  this  method  is  its 
applicability  to  analog  computers. 

TAEEUSZ  LESER  JOHN  T.  HARRISON 


\\ 
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FORAST  PROGRAM 


METHOD  OF  LIMES 

PROGRAMMER-  J.  T.  HARRISON 


COMM  GIVEN- 


COMM  DEL  U  *  0  IN  R 

COMM  U(X,0)  *  U(X,6)  *  lOOXf  12-X ) 

COMM  U(0,Y)  *  U(12iY)  *  100Y(8-Y) 

COMM  FIND  U(X«Y)  IN  REGION  R 

COMM  NOTATION- 


COMM  H  *  STEP  SIZE 

COMM  N  *  NUMBER  OF  LINES  TO  FINO  (0<Y<8). 

COMM  L  *  LENGTH  OF  Y  (0<*Y<*8). 

COMM  B1  *  Ul(O)  *  1200  -  B4  *  Ul( 12)  *  1200 

COMM  B2  *  U2(0)  *  1600  •  B5  *  U2I12)  *  1600 

COMM  B3  *  U3(0)  *  1200  9  86  *  U3<12)  *  1200 

COMM  Y  *  X 

COMM  YI  *  UK(X)  .  ( 1*1,2, 3)  *  (K*l,2,3). 

COMM  YI  *  VK( X)  *  U«K(X),(I*4,5,6)  »  (K*l,2,3). 

COMM  Y'l  *  0(UK)/DX  »  (1*1, 2, 3)  ,  (K*l,2,3). 

COMM  Y • I  *  D( VK)/DX  ,  (1*4,5,61  ,  (K-1,2,3). 


COMM 

COMM 

COMM 

COMM 

COMM 


Ml , 1  -  2N  X  2N  MATRIX  TO  FINE  THE  C*S. 

Cl  -  (1*1,2, ...6)  CONSTANTS  TO  BE  DETERMINED. 
CUI  -  ( 1*1,2,. ..21)  FINAL  SOLUTIONS  FOR  UK(X). 
Ql  -  (1*0,1,. ..6)  ERROR  TERMS  FOR  SUBROUTINE. 
UI  -  (1*1,2,... 230)  TEMP.  STORAGE. 
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BL0C!Y-V6HY*~Y*6MQ-Q6)(B1-B6)S 
BLOC  C  U1-U230 ) ( CU 1-CU45 M  Ml • 1-M6  «  7 ) S 
BLOCCC1-C6IS 
SVN  !X*Y|S 
DELX  DEC  C.m 

START  INHH-2IS  INT(N*3)S  I  NT  ( L-8 ) « 

PRINT-FORMAT  I F3  )**<  H  *  XHK  N  *  ><N) 

CONT<  L  *  XUS 

ENTER(PRINTB)S  ENTER! PRINTS )S 

COMM  BOUNDARY  CONDITIONS 

B1*1200S  B2*1600S  83*12001 
B4*1200S  85*16001  B6«1200S 

COMM  INTERGRATE  !0<*X<*12)»  BY  MEANS  OF 

COMM  A  SUBROUTINE#  RUNGE-KUTTA  GILL  TO  APPROXIMATE 

COMM  0R01NARY  01FF.  EQS.  FROM  7  INITIAL  CONDITIONS 

Y*«1S 

EPS*0ELX*.5S 
SET! TC*0) !  1*0)1 

1.0  READ-FORMAT ( FI)-! 7) NOS .AT! YO IS  INC!TC*TC«1)S 

SET!C*0)S 

U1#I*YIS  U2«I*Y2S  U3.I-Y3S 
INC! I*I+3)S 

ENTER! R.K.G.  )  I0ELX )  (7MEVAL*  Y)  ( Y)  ! Y' ) !Q)S 
COUNT ! 20  ) IN(C)G0T0!R.K.G1)S  G0T0(4.0)S 

COMM  EVALUATE  THE  V*S. 

EVAL* Y  Y* 1*Y4S  Y»2*Y5S  Y*3*Y6S 

IF- I NT! TC*7)G0T0!2.0)S 
HOM  *  OS  GOTO!  3.0IS 
2.0  HOM  *  20-30*X( 12-XIS 

3.0  AA*.6*V1-.3«Y24H0MS  BB*-.3»!  Y3-2*Y2m  )S 

CC*.6»Y3-.3*Y2^H0MS 
V*l*(99»AA-10*BB*CC)/98S 
V*2*« 100*BB-10*CC-10*AAI/98S 
V*3*!AA-10«BB»99«CC)/98S 
Y*4*V* IS  Y'5*V,2S  Y*6*V*3S 

GOTOIR.K.GDIS 
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COMM 


STORE  DATA 


4.0 


6.0 


FINDC 

7.0 


8.0 

8.1 

8.2 

8.3 

9.0 

10.0 


11.0 

12.0 


FI 

F2 

F3 


U1«I-Y1S  U2« I*Y2X  U3» I*Y3X 
INC(IM*3)X 

IF! X* 12) WITHIN I  EPS I GOTO! 6.0)X 

SET(C»0)X  G0T0(R.K.G1 )X 

IF-INTI TC*7)G0T0(FIN0C}X  GOTO( 1.0 )X 


COMM  DETERMINE  THE  C*S  FROM  UK(O)  AND  UKI12). 
COMM  FORM  A  2N  X  2N  MATRIX. 


SET! 1 1*0) ( JJ-O) (T*0) I J*0) (KK*0)X 
Ml* 1* I l*Ul * JJX  INC(JJ*JJ+21)(II*II+1)X 
C0UNTI6) 1N(T)G0T0(7.0)X 

Ml  *  1  *  I 1*81 *KK-U1 * JJX  INC(KK*KM1)  ( 1I*II«1)S 
SET ( T*0)X 

COUNT (3) INI J) GOTO (8.0 IS  GOTO (9.0) I 
IF-INT(KK<*2)G0T0(8.1)1  G0T0I8.2IX 
SET ( JJ*0)X  G0T0I8.3JX 
SET ( J J*18) X 

INT ( J J* J J+J ) X  GOTO! 7.0)X 

IF-INT(KK*6)G0T0(10.0)X  SET(J*0)X  G0T0I8.0IX 
ENTER! S.N.E. ) (M1*1)(6)(C1)X 


COMM  FINO  UK( X )  -|K*1,2.3)  IN  THE  REGION  R. 


SET ( K*0 ) ( 1 1*0) ( J  J*0 ) 1 1*0) IP* 126) X 

CLEAR (45) NOS. ATI CU1 )X  SET ( KK*1 ) ( CT*0 ) X 

PRINTOl  1XY*2>9XY*4>9XY*6>X 

ENTER! PRINTB ) X 

INT (KI*K+I ) X  SET! 1 1*0) X 

CU1*KI«C1* I I*U1* JJ+CU1*KIX 

I  NCI JJ>JJ«21)X 

COUNT (6) INI  1 1 ) GOTO  1 12.0) X 

CU1*KI*CU1 *KI *U1 *PX  INCI P*P*i )X 

INTI JJ«I^KK)X 


COUNT (3)IN(I)G0T0(11.0)X 
PRINT-FORMAT!  F2)-<X*XCT)  I  3 ) NOS.  ATICUl.X )X 
SET  1 1  *0)  X  INC<K«KO)(CT*CT*2)!KK»IMU3)X 
IF- 1  NT  I  K*21 )  GOTO  IN.  PROS )  X  ENTER!  PRINTS  )X  G0T0(U.0)X 
GOTO! 11.0) X 
FORM! 10-10)1 1-7 )X 

FORM 1 4-3) 1 3-2) (1-1)1 12-4-10) 1 3-2 ) 1 1-3 )X 
FORM I 4-3) I 1-3)X 
ENO  GOTO I  START)! 
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NETHOD  OF  LINES 


INPUT 


Y*X 

VI 

Y2 

Y3 

Y4 

Y5 

Y6 

O. 

u 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

1. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

1. 

0* 

0. 

0. 

0. 

0* 

0. 

0. 

1* 

0. 

0. 

0. 

0. 

0* 

0. 

0* 

I. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

1. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

METHOO  OF  LINES 


OUTPUT 

H  *  2  N  *  3  L  *  8 


Y»2 

V*4 

V*6 

X* 

1200.0000 

1600.0000 

1200.0000 

X* 

2 

1907.5268 

1944.4311 

1907.5268 

X* 

4 

2581.9081 

2420.4529 

2581.9081 

x= 

6 

2859.0200 

2620.3148 

2839.0200 

X* 

8 

2581.9081 

2420.4529 

2581.9081 

x» 

10 

1907.5268 

1944.4311 

1907.5268 

X* 

12 

1200.0000 

1600.0000 

1200.0000 
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